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Abstract 

In order to account for the observable Universe, any comprehensive 
theory or model of cosmology must draw from many disciplines of physics, 
including gauge theories of strong and weak interactions, the hydrody- 
namics and microphysics of baryonic matter, electromagnetic fields, and 
spacetime curvature, for example. Although it is difficult to incorporate 
all these physical elements into a single complete model of our Universe, 
advances in computing methods and technologies have contributed signif- 
icantly towards our understanding of cosmological models, the Universe, 
and astrophysical processes within them. A sample of numerical calcula- 
tions (and numerical methods) applied to specific issues in cosmology are 
reviewed in this article: from the Big Bang singularity dynamics to the 
fundamental interactions of gravitational waves; from the quark-hadron 
phase transition to the large scale structure of the Universe. The em- 
phasis, although not exclusively, is on those calculations designed to test 
different models of cosmology against the observed Universe. 
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1 Introduction 



Numerical investigations of cosmological spacetimes can be categorized into 
two broad classes of calculations, distinguished by their computational (or even 
philosophical) goals: 1) geometrical and mathematical principles of cosmolog- 
ical models, and 2) physical and astrophysical cosmology. In the former, the 
emphasis is on the geometric framework in which astrophysical processes occur, 
namely the cosmological expansion, shear, and singularities of the many models 
allowed by the theory of general relativity. In the latter, the emphasis is on the 
cosmological and astrophysical processes in the real or observable Universe, and 
the quest to determine the model which best describes our Universe. The former 
is pure in the sense that it concerns the fundamental nonlinear behavior of the 
Einstein equations and the gravitational field. The latter is more complex as 
it addresses the composition, organization, and dynamics of the Universe from 
the small scales (fundamental particles and elements) to the large (galaxies and 
clusters of galaxies). However the distinction is not always so clear, and geo- 
metric effects in the spacetime curvature can have significant consequences for 
the evolution and observation of matter distributions. 

Any comprehensive model of cosmology must therefore include nonlinear 
interactions between different matter sources and spacetime curvature. A real- 
istic model of the Universe must also cover large dynamical spatial and tempo- 
ral scales, extreme temperature and density distributions, and highly dynamic 
atomic and molecular matter compositions. In addition, due to all the var- 
ied physical processes of cosmological significance, one must draw from many 
disciplines of physics to model curvature anisotropies, gravitational waves, elec- 
tromagnetic fields, nucleosynthesis, particle physics, hydrodynamic fluids, etc. 
These phenomena are described in terms of coupled nonlinear partial differen- 
tial equations and must be solved numerically for general inhomogeneous space- 
times. The situation appears extremely complex, even with current technologi- 
cal and computational advances. As a result, the codes and numerical methods 
that have been developed to date are designed to investigate very specific prob- 
lems with either idealized symmetries or simplifying assumptions regarding the 
metric behavior, the matter distribution/composition or the interactions among 
the matter types and spacetime curvature. 

It is the purpose of this article to review published numerical cosmological 
calculations addressing problems from the very early Universe to the present; 
from the purely geometrical dynamics of the initial singularity to the large 
scale structure of the Universe. There are three major sections: § || where a 
brief overview is presented of various defining events occurring throughout the 
history of our Universe and in the context of the standard model; § || where 
summaries of early Universe and relativistic cosmological calculations are pre- 
sented; and § ^ which focuses on structure formation in the post-recombination 
epoch and on testing cosmological models against observations. Following a 
few conclusion statements in § H, an appendix § ra discusses the basic Einstein 
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equations, kinematic considerations, matter source equations with curvature, 
and the equations of pcrturbative physical cosmology on background isotropic 
models. References to numerical methods are also supplied and reviewed for 
each case. 
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2 Background 



2.1 A brief chronology 

With current observational constraints, the physical state of our Universe, as 
understood in the context of the standard or Friedmann-Lemaitre-Robertson- 
Walker (FLRW) model, can be crudely extrapolated back to ~ 1CP 34 seconds 
after the Big Bang, before which the classical description of general relativity 
is expected to give way to a quantum theory of gravity. At the earliest times, 
the Universe was a plasma of relativistic particles consisting of quarks, leptons, 
gauge bosons, and Higgs bosons represented by scalar fields with interaction 
and symmetry regulating potentials. It is believed that several spontaneous 
symmetry breaking (SSB) phase transitions occured in the early Universe as 
it expanded and cooled, including the grand unification transition (GUT) at 
~ 1CU 34 seconds after the Big Bang in which the strong nuclear force split off 
from the weak and electromagnetic forces (this also marks an era of inflation- 
ary expansion and the origin of matter-antimatter asymmetry through baryon, 
charge conjugation, and charge + parity violating interactions and nonequi- 
librium effects); the electroweak (EW) SSB transition at ~ 1CP 11 s when the 
weak nuclear force split from the electromagnetic force; and the chiral or quan- 
tum chromodynamic (QCD) symmetry breaking transition at ~ 10 -5 s during 
which quarks condensed into hadrons. The most stable hadrons (baryons, or 
protons and neutrons comprised of three quarks) survived the subsequent period 
of baryon-antibaryon annihilations, which continued until the Universe cooled 
to the point at which new baryon-antibaryon pairs could no longer be produced. 
This resulted in a large number of photons and relatively few surviving baryons. 
A period of primordial nucleosynthesis followed from ~ 1CU 2 to ~ 10 2 s during 
which light element abundances were synthesized to form 24% helium with trace 
amounts of deuterium, tritium, helium-3, and lithium. 

By ~ 10 11 s, the matter density became equal to the radiation density as 
the Universe continued to expand, identifying the start of the current matter- 
dominated era and the beginning of structure formation. Later, at ~ 10 13 s 
(3 x 10 5 years), the free ions and electrons combined to form atoms, effectively 
decoupling the matter from the radiation field as the Universe cooled. This de- 
coupling or post-recombination epoch marks the surface of last scattering and 
the boundary of the observable (via photons) Universe. Assuming a hierar- 
chical Cold Dark Matter (CDM) structure formation scenario, the subsequent 
development of our Universe is characterized by the growth of structures with 
increasing size. For example, the first stars are likely to have formed at t ~ 10 8 
years from molecular gas clouds when the Jeans mass of the background bary- 
onic fluid was approximately 10 4 M a , as indicated in Figure ^. This epoch of 
pop III star generation is followed by the formation of galaxies at t ~ 10 9 years 
and subsequently galaxy clusters. Though somewhat controversial, estimates of 
the current age of our Universe range from 10 to 20 Gyrs, with a present-day 
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linear structure scale radius of about 8 h^ 1 Megaparsecs, where h is the Hubble 
parameter (compared to 2-3 Megaparsecs typical for the virial radius of rich 
galaxy clusters). 

Figure 1: Schematic depicting the general sequence of events in the post- 
recombination Universe. The solid and dotted lines potentially track the Jeans 
mass of the average baryonic gas component from the recombination epoch at 
z ~ 10 3 to the current time. A residual ionization fraction o/nH+/ n H ~ 10~ 4 
following recombination allows for Compton interactions with photons to z ~ 
200, during which the Jeans mass remains constant at 1O 5 M0. The Jeans mass 
then decreases as the Universe expands adiabatically until the first collapsed 
structures form sufficient amounts of hydrogen molecules to trigger a cooling 
instability and produce pop III stars at z ~ 20. Star formation activity can 
then reheat the Universe and raise the mean Jeans mass to above 10° Mq. This 
reheating could affect the subsequent development of structures such as galaxies 
and the observed Lyman-alpha clouds. 



2.2 Successes of the standard model 

The isotropic and homogeneous FLRW cosmological model has been so suc- 
cessful in describing the observable Universe that it is commonly referred to as 
the "standard model". Furthermore, and to its credit, the model is relatively 
simple so that it allows for calculations and predictions to be made of the very 
early Universe, including primordial nucleosynthesis at 10~ 2 seconds after the 
Big Bang, and even particle interactions approaching the Planck scale at 10~ 43 
seconds. At present, observational support for the standard model includes: 

• the expansion of the Universe as verified by the redshifts in galaxy spectra 
and quantified by measurements of the Hubble constant Hq = lOO/i km s _1 
Mpc" 1 ; 

• the deceleration parameter observed in distant galaxy spectra (although 
uncertainties about galactic evolution, intrinsic luminosities, and standard 
candles prevent an accurate estimate); 

• the large scale isotropy and homogeneity of the Universe based on temper- 
ature anisotropy measurements of the microwave background radiation 
and peculiar velocity fields of galaxies (although the light distribution 
from bright galaxies is somewhat contradictory); 

• the age of the Universe which yields roughly consistent estimates between 
the look-back time to the Big Bang in the FLRW model and observed 
data such as the oldest stars, radioactive elements, and cooling of white 
dwarf stars; 
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• the cosmic microwave background radiation suggests that the Universe 
began from a hot Big Bang and the data is consistent with a mostly 
isotropic model and a black body at temperature 2.7 K; 

• the abundance of light elements such as 2 H, 3 He, 4 He and 7 Li, as predicted 
from the FLRW model, is consistent with observations and provides a 
bound on the baryon density and baryon-to-photon ratio; 

• the present mass density, as determined from measurements of luminous 
matter and galactic rotation curves, can be accounted for by the FLRW 
model with a single density parameter (f2o) to specify the metric topology; 

• the distribution of galaxies and larger scale structures can be reproduced 
by numerical simulations in the context of inhomogeneous perturbations 
of the FLRW models. 



Because of these successes, most work in the field of physical cosmology (see 
§ ||) has utilized the standard model as the background spacetime in which the 
large scale structure evolves, with the ambition to further constrain parameters 
and structure formation scenarios through numerical simulations. The reader is 
referred to [§4j for a more in-depth review of the standard model, and to [102, 



119 1 for a summary of observed cosmological parameter constraints and best fit 



"concordance" models. 
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3 Relativistic Cosmology 



This section is organized to track the chronological events in the history of the 
early or relativistic Universe, focusing mainly on four defining moments: 1) the 
Big Bang singularity and the dynamics of the very early Universe; 2) inflation 
and its generic nature; 3) QCD phase transitions; and 4) primordial nucleosyn- 
thesis and the freeze-out of the light elements. The late or post-recombination 
epoch is reserved to a separate section § f|. 



3.1 Singularities 

3.1.1 Mixmaster dynamics 

Belinsky, Lifshitz and Khalatnikov (BLK) ||(| and Misner p5j discovered 
that the Einstein equations in the vacuum homogeneous Bianchi type IX (or 
Mixmaster) cosmology exhibit complex behavior and are sensitive to initial con- 
ditions as the Big Bang singularity is approached. In particular, the solutions 
near the singularity are described qualitatively by a discrete map [ p7| p9| rep- 
resenting different sequences of Kasner spacetimes 

ds 2 = -dt 2 + t 2pi dx 2 + t 2p2 dy 2 + t 2p3 dz 2 , (1) 

with time changing exponents pi , but otherwise constrained by pi + P2 + P3 = 
Pi + P% +P§ = 1- Because this discrete mapping of Kasner epochs is chaotic, 
the Mixmaster dynamics is presumed to be chaotic as well. 



Figure 2: Contour plot of the Bianchi type IX potential V , where f}± are the 
anisotropy canonical coordinates. Seven level surfaces are shown at equally 
spaced decades ranging from 1CU 1 to 10 5 . For large isocontours (V > 1 ), the 
potential is open and exhibits a strong triangular symmetry with three narrow 
channels extending to spatial infinity. For V < 1, the potential closes and is 
approximately circular for f3± <C 1 . 



Mixmaster behavior can be studied in the context of Hamiltonian dynamics, 
with a Hamiltonian [M 



in 



-pft+pi+pt+e^iV-l), 



(2) 



and a semi-bounded potential arising from the spatial scalar curvature (whose 
level curves are plotted in Figure 0) 



V = l + -e-^+ + -e 4 ^ 



cosh(4V3/3_) - 1 - -e- 2f3 + cosh(2V3/3-), (3) 
o 



where e" and (3± are the scale factor and anisotropics, and p a and p± are the 
corresponding conjugate variables. A solution of this Hamiltonian system is 
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an infinite sequence of Kasner epochs with parameters that change when the 
phase space trajectories bounce off the potential walls, which become exponen- 
tially steep as the system evolves towards the singularity. Several numerical 
calculations of the dynamical equations arising from (^|) and (^|) have indicated 
that the Liapunov exponents of the system vanish, in apparent contradiction 
with the discrete maps ff^ , and putting into question the characterization 
of Mixmaster dynamics as chaotic. However, it has since been shown that the 
usual definition of the Liapunov exponents is ambiguous in this case as it is not 
invariant under time reparametrizations, and that with a different time variable 
one obtains positive exponents [^3[ [35). Also, coordinate independent methods 
using fractal basin boundaries to map basins of attraction in the space of initial 
conditions indicates Mixmaster spacetimes to be chaotic . 

Although BLK conjectured that local Mixmaster oscillations might be the 
generic behavior for singularities in more general classes of spacetimes |3Cj, it 
is only recently that this conjecture has begun to be supported by numerical 



evidence (see Section 3.1.2 and [pi 



3.1.2 AVTD vs. BLK oscillatory behavior 

As noted in § |3.1.l[ an interesting and important issue in classical cosmology 
is whether or not the generic Big Bang singularity is locally of a Mixmaster or 
BLK type, with complex oscillatory behavior as the singularity is approached. 
Most of the Bianchi models, including the Kasner solutions (0), are charac 



terized by either open or no potentials in the Hamiltonian framework [ 109 1 , 
and exhibit essentially monotonic or Asymptotically Velocity Term Dominated 
(AVTD) behavior. 

Considering inhomogeneous spacetimes, Isenberg and Moncrief ]8l| proved 
that the singularity in the polarized Gowdy model is AVTD type, as are more 
general polarized T 2 symmetric cosmologies [ p4[ . Early numerical studies using 
symplectic methods have confirmed these conjectures and found no evidence of 
BLK oscillations, even in T 3 x R spacetimes with U(l) symmetry j32| (although 
there were concerns about the solutions due to difficulties in resolving steep 
spatial gradients near the singularity Q), which were addressed later by Hern 



and Stewart [jf5j for the Gowdy T 3 models). However, Weaver et al. \124\ have 
established the first evidence through numerical simulations that Mixmaster dy- 
namics can occur in (at least a restricted class of) inhomogeneous spacetimes 
which generalize the Bianchi type VIo with a magnetic field and two-torus sym- 
metry. More recently, Berger and Moncrief (3^, [57j have shown U(l) symmet- 
ric vacuum cosmologies to exhibit local Mixmaster dynamics, which tends to 
support the BLK conjecture. Despite numerical difficulties in resolving steep 
gradients (which they managed by enforcing the Hamiltonian constraint and 
spatially averaging the solutions), Berger and Moncrief have confirmed their 
findings under increased spatial resolution and changes in initial data. 
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3.2 Inflation 



The inflation paradigm is frequently invoked to explain the flatness (flo w 1 in 
the context of the FLRW model) and nearly isotropic nature of the Universe 
at large scales, attributing them to an era of exponential expansion at about 
10~ 34 s after the Big Bang. This expansion acts as a strong dampening mecha- 
nism to random curvature or density fluctuations, and may be a generic attractor 
in the space of initial conditions. An essential component needed to trigger this 
inflationary phase is a scalar or inflaton field <f> representing spin zero particles. 
The vacuum energy of the field acts as an effective cosmological constant that 
regulates GUT symmetry breaking, particle creation, and the reheating of the 
Universe through an interaction potential V(<f>) derived from the form of sym- 
metry breaking that occurs as the temperature of the Universe decreases. Early 
analytic studies focused on simplified models, treating the interaction potential 
as flat near its local maximum where the field does not evolve significantly and 
where the formal analogy to an effective cosmological constant approximation 
is more precise. However, to properly account for the complexity of inflaton 



fields, the full dynamical equations (see § |6.3.2| ) must be considered together 
with consistent curvature, matter and other scalar field couplings. Also, many 
different theories of inflation and vacuum potentials have been proposed (see, 
for example, a recent review by Lyth and Riotto ]9l| and an introductory arti- 
cle by Liddle p(|), and it is not likely that simplified models can elucidate the 



full nonlinear complexity of scalar fields (see § 3.3) nor the generic nature of 
inflation. 

In order to study whether inflation can occur for arbitrary anisotropic and 
inhomogeneous data, many numerical simulations have been carried out with 
different symmetries, matter types and perturbations. A sample of such calcu- 
lations is described in the following paragraphs. 

3.2.1 Plane symmetry 

Kurki-Suonio ct al. |p5| ex tended the planar cosmology code of Centrella and 
Wilson |34|, |35| (see § |3.6|) to include a scalar field and simulate the onset of 
inflation in the early Universe with an inhomogeneous Higgs field and a perfect 
fluid with a radiation equation of state p = p/3, where p is the pressure and p 
is the energy density. Their results suggest that whether inflation occurs or not 
can be sensitive to the shape of the potential </>. In particular, if the shape is flat 
enough and satisfies the slow-roll conditions (essentially upper bounds on dV/ d(f> 
and d 2 V/d<f> 2 |84| near the false vacuum <f> ~ 0), even large initial fluctuations 
of the Higgs field do not prevent inflation. They considered two different forms 
of the potential: a Coleman- Weinberg type with interaction strength A 



V(4>) = A0 4 



.,2 
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(4) 



which is very flat close to the false vacuum and does inflate; and a rounder "</> 4 " 
type 

V{<P) = A(0 2 - a 2 ) 2 (5) 
which, for their parameter combinations, does not. 

3.2.2 Spherical symmetry 

Goldwirth and Piran J7l| studied the onset of inflation with inhomogeneous ini- 
tial conditions for closed, spherically symmetric spacetimes containing a massive 
scalar field and radiation field sources (described by a massless scalar field). In 
all the cases they considered, the radiation field damps quickly and only an 
inhomogeneous massive scalar field remains to inflate the Universe. They find 
that small inhomogeneities tend to reduce the amount of inflation and large ini- 
tial inhomogeneities can even suppress the onset of inflation. Their calculations 
indicate that the scalar field must have "suitable" initial values over a domain 
of several horizon lengths in order for inflation to begin. 



3.2.3 Bianchi V 

Anninos et al. |l2| investigated the simplest Bianchi model (type V) background 
that admits velocities or tilt in order to address the question of how the Universe 
can choose a uniform reference frame at the exit from inflation, since the de 
Sitter metric does not have a preferred frame. They find that inflation does 
not isotropize the Universe in the short wavelength limit. However, if inflation 
persists, the wave behavior eventually freezes in and all velocities go to zero 
at least as rapidly as tanh/3 ~ -R -1 , where f3 is the relativistic tilt angle (a 
measure of velocity), and R is a typical scale associated with the radius of the 
Universe. Their results indicate that the velocities eventually go to zero as 
inflation carries all spatial variations outside the horizon, and that the answer 
to the posed question is that memory is retained and the Universe is never really 
de Sitter. 



3.2.4 Gravity waves + cosmological constant 

In addition to the inflaton field, one can consider other sources of inhomogene- 
ity, such as gravitational waves. Although linear waves in de Sitter space will 
decay exponentially and disappear, it is unclear what will happen if strong 



waves exist. Shinkai & Maeda [116] investigated the cosmic no-hair conjecture 
with gravitational waves and a cosmological constant (A) in ID plane sym- 
metric vacuum spacetimes, setting up Gaussian pulse wave data with ampli- 
tudes 0.02A < max(VT) < 80A and widths 0.08 < I < 2.5 Z H , where I 
is the invariant constructed from the 3-Riemann tensor and Zh = y/S/A is 
the horizon scale. They also considered colliding plane waves with amplitudes 
40A < max(VI) < 125A and widths 0.08 Z H < I < 0.1 Z H - They find that for 
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any large amplitude or small width inhomogeneity in their parameter sets, the 
nonlinearity of gravity has little effect and the spacetime always evolves towards 
de Sitter. 



3.2.5 3D inhomogeneous spacetimes 

The previous paragraphs discussed results from highly symmetric spacetimes, 
but the possibility of inflation remains to be established for more general inho- 
mogeneous and nonperturbative data. In an effort to address this issue, Kurki- 
Suonio et al. p6| investigated fully three-dimensional inhomogeneous spacetimes 
with a chaotic inflationary potential V{<j>) = A0 4 /4. They considered basically 
two types of runs: small and large scale. In the small scale run, the grid length 
was initially set equal to the Hubble length so the inhomogeneities are well in- 
side the horizon and the dynamical time scale is shorter than the expansion or 
Hubble time. As a result, the perturbations oscillate and damp, while the field 
evolves and the spacetime inflates. In the large scale run, the inhomogeneities 
are outside the horizon and they do not oscillate. They maintain their shape 
without damping and, because larger values of 4> lead to faster expansion, the 
inhomogeneity in the expansion becomes steeper in time since the regions of 
large <fi and high inflation stay correlated. Both runs have sufficient inflation to 
solve the flatness problem. 



3.3 Chaotic scalar field dynamics 

Many studies of cosmological models generally reduce complex physical systems 
to simplified or even analytically integrable systems. In sufficiently simple mod- 
els the dynamical behavior (or fate) of the Universe can be predicted from a 
given set of initial conditions. However, the Universe is composed of many dif- 
ferent nonlinear interacting fields including the inflaton field which can exhibit 
complex chaotic behavior. For example, Cornish and Levin |57j consider a ho- 
mogenous isotropic closed FLRW model with various conformal and minimally 



coupled scalar fields (see § 6.3.2 ). They find that even these relatively simple 
models exhibit chaotic transients in their early pre-inflationary evolution. This 
behavior in exiting the Planck era is characterized by fractal basins of attraction, 
with attractor states being to (1) inflate forever, (2) inflate over a short period 
of time then collapse, or (3) collapse without inflating. Monerat et al. (9^] inves- 
tigated the dynamics of the pre-inflationary phase of the Universe and its exit to 
inflation in a closed FLRW model with radiation and a minimally coupled scalar 
field. They observe complex behavior associated with saddle-type critical points 
in phase space that give rise to deSitter attractors with multiple chaotic exits 
to inflation that depend on the structure of the scalar field potential. These 
results suggest that distinctions between exits to inflation may be manifested in 
the process of reheating and as a selected spectrum of inhomogeneous perturba- 
tions influenced by resonance mechanisms in curvature oscillations. This could 
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possibly lead to fractal patterns in the density spectrum, gravitational waves, 
CMBR field, or galaxy distribution that depend on specific details including the 
number of fields, the nature of the fields, and their interaction potentials. 

Chaotic behavior can also be found in more general applications of scalar 
field dynamics. Anninos et al. pl| investigated the nonlinear behavior of collid- 
ing kink-antikink solitons or domain walls described by a single real scalar field 
with self-interaction potential X(<j> 2 — l) 2 . Domain walls can form as topological 
defects during the spontaneous symmetry breaking period associated with phase 
transitions, and can seed density fluctuations in the large scale structure. For 
collisional time scales much smaller than the cosmological expansion, they find 
that whether a kink-antikink collision results in a bound state or a two-soliton 
solution depends on a fractal structure observed in the impact velocity param- 
eter space. The fractal structure arises from a resonance condition associated 
with energy exchanges between translational modes and internal shape-mode 
oscillations. The impact parameter space is a complex self-similar fractal com- 
posed of sequences of different n-bounce (the number of bounces or oscillations 
in the transient semi-coherent state) reflection windows separated by regions of 
oscillating bion states (see Figure ^) . They compute the Lyapunov exponents 
for parameters in which a bound state forms to demonstrate the chaotic nature 
of the bion oscillations. 

Figure 3: Fractal structure of the transition between reflected and captured states 
for colliding kink-antikink solitons in the parameter space of impact velocity 
for a \{(f> 2 — l) 2 scalar field potential. The top image (a) shows the 2-bounce 
windows in dark with the rightmost region (v/c> 0.25) representing the single- 
bounce regime above which no captured state exists, and the leftmost white region 
(v/c < 0.19) representing the captured state below which no reflection windows 
exist. Between these two marker velocities, there are 2-bounce reflection states 
of decreasing widths separated by regions of bion formation. Zooming in on the 
domain outlined by the dashed box, a self-similar structure is apparent in the 
middle image (b), where now the dark regions represent 3-bounce windows of 
decreasing widths. Zooming in once again on the boundaries of these 3-bounce 
windows, a similar structure is found as shown in the bottom image (c) but 
with ^-bounce reflection windows. This pattern of self- similarity with n-bounce 
windows is observed at all scales investigated numerically. 



3.4 Quark-hadron phase transition 

The standard picture of cosmology assumes that a phase transition (associ- 
ated with chiral symmetry breaking after the electroweak transition) occurred 
at approximately 10 -5 seconds after the Big Bang to convert a plasma of free 
quarks and gluons into hadrons. Although this transition can be of significant 
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cosmological importance, it is not known with certainty whether it is of first 
order or higher, and what the astrophysical consequences might be on the sub- 
sequent state of the Universe. For example, the transition may give rise to 
significant baryon number inhomogeneities which can influence the outcome of 
primordial nucleosynthesis as evidenced in the distribution and averaged light 
element abundances. The QCD transition and baryon inhomogeneities may also 
play a significant and potentially observable role in the generation of primordial 
magnetic fields. 



Rczolla et al. [106] considered a first order phase transition and the nucle- 
ation of hadronic bubbles in a supercooled quark-gluon plasma, solving the rel- 
ativistic Lagrangian equations for disconnected and evaporating quark regions 
during the final stages of the phase transition. They numerically investigated 
a single isolated quark drop with an initial radius large enough so that surface 
effects can be neglected. The droplet evolves as a self-similar solution until it 
evaporates to a sufficiently small radius that surface effects break the similarity 
solution and increase the evaporation rate. Their simulations indicate that, in 
neglecting long-range energy and momentum transfer (by electromagnetically 
interacting particles) and assuming that the baryon number is transported with 
the hydrodynamical flux, the baryon number concentration is similar to what 
is predicted by chemical equilibrium calculations. 



Kurki-Suonio and Lainc |87 studied the growth of bubbles and the decay of 
droplets using a spherically symmetric code that accounts for a phenomenolog- 
ical model of the microscopic entropy generated at the phase transition surface. 
Incorporating the small scale effects of the finite wall width and surface tension, 
but neglecting entropy and baryon flow through the droplet wall, they demon- 
strate the dynamics of nucleated bubble growth and quark droplet decay. They 
also find that evaporating droplets do not leave behind a global rarefaction wave 
to dissipate any previously generated baryon number inhomogeneity. 

3.5 Nucleosynthesis 

Observations of the light elements produced during Big Bang nucleosynthesis 
following the quark/hadron phase transition (roughly 10~ 2 -10 2 seconds after 
the Big Bang) are in good agreement with the standard model of our Universe 



(see § 2.2). However, it is interesting to investigate other more general models 



to assert the role of shear and curvature on the nucleosynthesis process. 



Rothman and Matzner [108 considered primordial nucleosynthesis in aniso- 
tropic cosmologies, solving the strong reaction equations leading to 4 He. They 
find that the concentration of 4 He increases with increasing shear due to time 
scale effects and the competition between dissipation and enhanced reaction 
rates from photon heating and neutrino blue shifts. Their results have been 
used to place a limit on anisotropy at the epoch of nucleosynthesis. Kurki- 
Suonio and Matzner [^8| extended this work to include 30 strong 2-particle 
reactions involving nuclei with mass numbers A < 7, and to demonstrate the 
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effects of anisotropy on the cosmologically significant isotopes 2 H, 3 He, 4 He and 
7 Li as a function of the baryon to photon ratio. They conclude that the effect 
of anisotropy on 2 H and 3 He is not significant, and the abundances of 4 He and 
7 Li increase with anisotropy in accord with |108| . 

Furthermore, it is possible that neutron diffusion, the process whereby neu- 
trons diffuse out from regions of very high baryon density just before nucle- 
osynthesis, can affect the neutron to proton ratio in such a way as to enhance 
deuterium and reduce 4 He compared to a homogeneous model. However, plane 
symmetric, general relativistic simulations with neutron diffusion [jS9| show that 
the neutrons diffuse back into the high density regions once nucleosynthesis be- 
gins there - thereby wiping out the effect. As a result, although inhomogeneities 
influence the element abundances, they do so at a much smaller degree then pre- 
viously speculated. The numerical simulations also demonstrate that, because of 
the back diffusion, a cosmological model with a critical baryon density cannot be 
made consistent with helium and deuterium observations, even with substantial 
baryon inhomogeneities and the anticipated neutron diffusion effect. 



3.6 Plane symmetric gravitational waves 

Gravitational waves are an inevitable product of the Einstein equations, and 
one can expect a wide spectrum of wave signals propagating throughout our 
Universe due to shear anisotropics, primordial metric and matter fluctuations, 
collapsing matter structures, ringing black holes, and colliding neutron stars, for 
example. The discussion here is restricted to the pure vacuum field dynamics 
and the fundamental nonlinear behavior of gravitational waves in numerically 
generated cosmological spacetimes. 

Centrella and Matzner |5^, studied a class of plane symmetric cosmologies 
representing gravitational inhomogeneities in the form of shocks or discontinu- 
ities separating two vacuum expanding Kasner cosmologies ([!]). By a suitable 
choice of parameters, the constraint equations can be satisfied at the initial time 
with an Euclidean 3-surface and an algebraic matching of parameters across the 
different Kasner regions that gives rise to a discontinuous extrinsic curvature 
tensor. They performed both numerical calculations and analytical estimates 
using a Green's function analysis to establish and verify (despite the numeri- 
cal difficulties in evolving discontinuous data) certain aspects of the solutions, 
including gravitational wave interactions, the formation of tails, and the singu- 
larity behavior of colliding waves in expanding vacuum cosmologies. 

Shortly thereafter, Centrella and Wilson |)4|, [5^] developed a polarized plane 
symmetric code for cosmology, adding also hydrodynamic sources with artificial 
visco sity m ethods for shock capturing and Barton's method for monotonic trans- 



port [126]. The evolutions are fully constrained (solving both the momentum 
and Hamiltonian constraints at each time step) and use the mean curvature slic- 
ing condition. This work was subsequently extended by Anninos et al. H [l0[ jof, 
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implementing more robust numerical methods, an improved parametric treat- 
ment of the initial value problem, and generic unpolarized metrics. 

In applications of these codes, Centrella |5l| investigated nonlinear gravity 
waves in Minkowski space and compared the full numerical solutions against 
a first order perturbation solution to benchmark certain numerical issues such 
as numerical damping and dispersion. A second order perturbation analysis 
was used to model the transition into the nonlinear regime. Anninos et al. Q 
considered small and large perturbations in the two degenerate Kasner models: 
Pi — P3 — or 2/3, and p2 = 1 or —1/3 respectively, where p>\ are parameters 
in the Kasner metric (]l|). Carrying out a second order perturbation expansion 
and computing the Newman-Penrose (NP) scalars, Riemann invariants and Bel- 
Robinson vector, they demonstrated, for their particular class of spacctimes, 
that the nonlinear behavior is in the Coulomb (or background) part represented 
by the leading order term in the NP scalar \&2, and not in the gravitational 
wave component. For standing-wave perturbations, the dominant second order 
effects in their variables are an enhanced monotonic increase in the background 
expansion rate, and the generation of oscillatory behavior in the background 
spacetime with frequencies equal to the harmonics of the first order standing- 
wave solution. Expanding their investigations of the Coulomb nonlinearity, 
Anninos and McKinney [Q used a gauge invariant perturbation formalism to 
construct constrained initial data for general relativistic cosmological sheets 
formed from the gravitational collapse of an ideal gas in a critically closed FLRW 



"background" model. Results are compared to the Newtonian Zel'dovich [ 128 1 
solution over a range of field strengths and flows. Also, the growth rates of 
nonlinear modes (in both the gas density and Riemann curvature invariants), 
their effect in the back-reaction to modify the cosmological scale factor, and 
their role in generating CMB anisotropies are discussed. 

3.7 Regge calculus model 

A unique approach to numerical cosmology (and numerical relativity in general) 
is the method of Regge Calculus in which spacetime is represented as a complex 
of 4-dimensional, geometrically flat simpliccs. The principles of Einstein's theory 
are applied directly to the simplicial geometry to form the curvature, action, and 
field equations, in contrast to the finite difference approach where the continuum 
field equations are differenced on a discrete mesh. 

A 3-dimensional code implementing Regge Calculus techniques was devel- 
oped recently by Gentle and Miller |5{| and applied to the Kasner cosmological 
model. They also describe a procedure to solve the constraint equations for time 
asymmetric initial data on two spacelike hypersurfaces constructed from tetra- 
hedra, since full 4-dimcnsional regions or lattices arc used. The new method 



is analogous to York's procedure (see [127] and § 3.4 ) where the conformal 
metric, trace of the extrinsic curvature, and momentum variables are all freely 
specifiable. These early results are promising in that they have reproduced the 



15 



continuum Kasner solution, achieved second order convergence, and sustained 
numerical stability. Also, Barnett et al. [^6| discuss an implicit evolution scheme 
that allows local (vertex) calculations for efficient parallelism. However, the 
Regge Calculus approach remains to be developed and applied to more general 
spacetimes with complex topologies, extended degrees of freedom, and general 
source terms. 
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4 Physical Cosmology 



The phrase "physical cosmology" is generally associated with the large (galaxy 
and cluster) scale structure of the post-recombination epoch where gravitational 
effects are modeled approximately by Newtonian physics on a uniformly ex- 
panding, matter dominated FLRW background. A discussion of the large scale 
structure is included in this review since any viable model of our Universe which 
allows a regime where strongly general relativistic effects are important must 
match onto the weakly relativistic (or Newtonian) regime. Also, since certain 
aspects of this regime are directly observable, one can hope to constrain or rule 
out various cosmological models and/or parameters, including the density (f2o), 
Hubble (Ho — 100 h km s _1 Mpc -1 ), and cosmological (A) constants. 

Due to the vast body of literature on numerical simulations of the post- 
recombination epoch, it is possible to mention only a very small fraction of all 
the published papers. Hence, the following summary is limited to cover just a 
few aspects of computational physical cosmology, and in particular those that 
can potentially be used to discriminate between cosmological model parameters, 
even within the realm of the standard model. 



4.1 Cosmic microwave background 

The Cosmic Microwave Background Radiation (CMBR), which is a direct relic 
of the early Universe, currently provides the deepest probe of cosmological struc- 
tures and imposes severe constraints on the various proposed matter evolution 
scenarios and cosmological parameters. Although the CMBR is a unique and 
deep probe of both the thermal history of the early Universe and the primor- 
dial perturbations in the matter distribution, the associated anisotropics are 
not exclusively primordial in nature. Important modifications to the CMBR 
spectrum can arise from large scale coherent structures, even well after the pho- 
tons decouple from the matter at redshift z ~ 10 3 , due to the gravitational 
redshifting of the photons through the Sachs- Wolfe effect arising from potential 



gradients 1 1 1 1 , O 



_ = $ e - $ r - / dt, (6) 

T J r a 

where the integral is evaluated from the emission (e) to reception (r) points along 
the spatial photon paths I, <£> is the gravitational potential, AT/T defines the 
temperature fluctuations, and a(t) is the cosmological scale factor in the stan- 
dard FLRW metric. Also, if the intergalactic medium (IGM) reionizes sometime 
after the decoupling, say from an early generation of stars, the increased rate of 
Thomson scattering off the free electrons will erase sub-horizon scale tempera- 
ture anisotropics, while creating secondary Doppler shift anisotropics. To make 
meaningful comparisons between numerical models and observed data, these 



effects (and others, see for example § 4.1.3 and references |79|, |82fl) must be in- 
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corporated self-consistently into the numerical models and to high accuracy in 
order to resolve the weak signals. 



4.1.1 Ray-tracing 

Many computational analyses based on linear perturbation theory have been 
carried out to estimate the temperature anisotropics in the sky (for example 
see |}2j and the references cited in |7£j). Although such linearized approaches 
yield reasonable results, they are not well-suited to discussing the expected 
imaging of the developing nonlinear structures in the microwave background. 
An alternative ray-tracing approach has been developed by Anninos et al. [jl3| 
to introduce and propagate individual photons through the evolving nonlinear 
matter structures. They solve the geodesic equations of motion and subject the 
photons to Thomson scattering in a probabilistic way and at a rate determined 
by the local density of free electrons in the model. Since the temperature fluc- 
tuations remain small, the equations of motion for the photons are treated as 
in the linearized limit, and the anisotropies are computed according to 

— = — (7) 

where 

1 + * = S?% (8) 

(fc^)r 

and the photon wave vector and matter rest frame four-velocity it M are eval- 
uated at the emission (e) and reception (r) points. Applying their procedure 
to a Hot Dark Matter (HDM) model of structure formation, Anninos et al. jl3| 
find the parameters for this model are severely constrained by COBE data such 
that floh 2 rs 1, where fio and h are the density and Hubble parameters. 

4.1.2 Effects of reionization 

In models where the IGM does not reionize, the probability of scattering after 
the photon-matter decoupling epoch is low, and the Sachs- Wolfe effect domi- 
nates the anisotropies at angular scales larger than a few degrees. However, if 
reionization occurs, the scattering probability increases substantially and the 
matter structures, which develop large bulk motions relative to the comoving 
background, induce Doppler shifts on the scattered CMBR photons and leave an 
imprint of the surface of last scattering. The induced fluctuations on subhori- 
zon scales in reionization scenarios can be a significant fraction of the primordial 



anisotropies, as observed by Tuluie et al. [ 1 22 1 . They considered two possible 
scenarios of reionization: A model that suffers early and gradual (EG) reion- 
ization of the IGM as caused by the photoionizing UV radiation emitted by 
decaying neutrinos, and the late and sudden (LS) scenario as might be applica- 
ble to the case of an early generation of star formation activity at high redshifts. 
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Considering the HDM model with Oo = 1 and h = 0.55, which produces CMBR 
anisotropics above current COBE limits when no reionization is included (see 
§ [4.1. lty , they find that the EG scenario effectively reduces the anisotropics to 
the levels observed by COBE and generates smaller Doppler shift anisotropies 
than the LS model, as demonstrated in Figure ||. The LS scenario of reioniza- 
tion is not able to reduce the anisotropy levels below the COBE limits, and can 
even give rise to greater Doppler shifts than expected at decoupling. 



Figure 4: The top two images represent temperature fluctuations (i.e., AT/T) 
due to the Sachs- Wolfe effect and Doppler shifts in a standard critically closed 
HDM model with no reionization and baryon fractions 0.02 (plate 1: 4° x 4°, 
rms = 2.8 x lCT 5 j and 0.2 (plate 2: 8° x 8°, rms = 3.4 x lCT 5 j. The bottom 
two plates image fluctuations in an "early and gradual" reionization scenario of 
decaying neutrinos with baryon fraction 0.02 (plate 3: 4° x4°, rms = 1.3 x 10~ 5 ; 
and plate 4: 8° x 8°, rms = 1.4 x 1(T 5 J. 



4.1.3 Secondary anisotropies 

Additional sources of CMBR anisotropy can arise from the interactions of pho- 
tons with dynamically ev olvin g matter structures and nonstatic gravitational 



potentials. Tuluie et al. [121| considered the impact of nonlinear matter con- 
densations on the CMBR in f2o < 1 Cold Dark Matter (CDM) models, focusing 
on the relative importance of secondary temperature anisotropies due to three 
different effects: 1) time-dependent variations in the gravitational potential of 
nonlinear structures as a result of collapse or expansion (the Rees-Sciama ef- 
fect); 2) proper motion of nonlinear structures such as clusters and superclusters 
across the sky; and 3) the decaying gravitational potential effect from the evo- 
lution of perturbations in open models. They applied the ray-tracing procedure 
of |H| to explore the relative importance of these secondary anisotropies as 
a function of the density parameter and the scale of matter distributions. 
They find that secondary temperature anisotropies are dominated by the de- 
caying potential effect at large scales, but that all three sources of anisotropy 
can produce signatures of order AT/T ~ 10 -6 as shown in Figure 0. 

In addition to the effects discussed here, many other sources of secondary 
anisotropies (such as gravitational lensing, the Vishniac effect accounting for 
matter velocities and flows into local potential wells, and the Sunyaev-Zel'dovich 



(§ 4.5.4 ) distortions from the Compton scattering of CMB photons by electrons 
in the hot cluster medium) can also be significant. See reference (7j| for a 
more complete list and thorough discussion of the different sources of CMBR 
anisotropies. 
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Figure 5: The top two images represent the proper motion and Rees-Sciama 
effects in the CMBR for a critically closed CDM model (upper left), together 
with the corresponding column density of voids and clusters over the same re- 
gion (upper right). The bottom two images show the secondary anisotropics 
dominated here by the decaying potential effect in an open cosmological model 
(bottom left), together with the corresponding gravitational potential over the 
same region (bottom right). The rms fluctuations in both cases are on the order 
o/±5 x 10~ 7 , though the open model carries a somewhat larger signature. 

4.2 Gravitational lensing 

Observations of gravitational lenses |112| provide measures of the abundance and 
strength of nonlinear potential fluctuations along the lines of sight to distant 
objects. Since these calculations are sensitive to the gravitational potential, they 
may be more reliable than galaxy and velocity field measurements as they are 
not subject to the same ambiguities associated with biasing effects. Also, since 
different cosmological models predict different mass distributions, especially at 
the higher redshifts, lensing calculations can potentially be used to confirm or 
discard competing cosmological models. 

As an alternative to solving the computationally demanding lens equations, 
Cen et al. [[f9) developed an efficient scheme to identify regions with surface 
densities capable of generating multiple images accurately for splittings larger 
than three arcseconds. They applied this technique to a standard CDM model 
with Oo = 1, and found that this model predicts more large angle splittings (> 
8") than are known to exist in the observed Universe. Their results suggest that 
the standard CDM model should be excluded as a viable model of our Universe. 
A similar analysis for a flat low density CDM model with a cosmological constant 
(Qo = 0.3, A/3-ffo = 0.7) suggests a lower and perhaps acceptable number of 
lensing events. However, an uncertainty in their studies is the nature of the 
lenses at and below the resolution of the numerical grid. They model the lensing 
structures as simplified Singular Isothermal Spheres (SIS) with no distinctive 
cores. 

Large angle splittings are generally attributed to larger structures such as 
clusters of galaxies, and there are indications that clusters have small but finite 
core radii r corc ~ 20 — 30 h^ 1 kpc. Core radii of this size can have an important 
effect on the probability of multiple imaging. Flores and Primack [ [S6| consid- 
ered the effects of assuming two different kinds of splitting sources: isothermal 
spheres with small but finite core radii p oc (r 2 + r 2 ^)^ 1 , and spheres with a 
Hcrnquist density profile p cx r _1 (r + a)~ 3 , where r corc ~ 20 — 30 h~ x kpc and 
a r-j 300 h~ x kpc. They find that the computed frequency of large-angle split- 
tings, when using the nonsingular profiles, can potentially decrease by more than 
an order of magnitude relative to the SIS case and can bring the standard CDM 
model into better agreement with observations. They stress that lensing events 
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are sensitive to both the cosmological model (essentially the number density of 
lenses) and to the inner lens structure, making it difficult to probe the models 
until the structure of the lenses, both observationally and numerically, is better 
known. 



4.3 First star formation 

In CDM cosmogonies, the fluctuation spectrum at small wavelengths has a loga- 
rithmic dependence at mass scales smaller than 10 s solar masses, which indicates 
that all small scale fluctuations in this model collapse nearly simultaneously in 
time. This leads to very complex dynamics during the formation of these first 
structures. Furthermore, the cooling in these fluctuations is dominated by the 
rotational/ vibrational modes of hydrogen molecules that were able to form us- 
ing the free electrons left over from recombination and those produced by strong 
shock waves as catalysts. The first structures to collapse may be capable of pro- 
ducing pop III stars and have a substantial influence on the subsequent thermal 
evolution of the intergalactic medium, as suggested by Figure [j], due to the ra- 
diation emitted by the first generation stars as well as supernova driven winds. 
To know the subsequent fate of the Universe and which structures will survive 
or be destroyed by the UV background, it is first necessary to know when and 
how the first stars formed. 



Ostriker and Gnedin [ 101 have carried out high resolution numerical simu- 
lations of the reheating and reionization of the Universe due to star formation 
bursts triggered by molecular hydrogen cooling. Accounting for the chemistry of 
the primeval hydrogen/helium plasma, self-shielding of the gas, radiative cool- 
ing, and a phenomenological model of star formation, they find that two distinct 
star populations form: the first generation pop III from Hi cooling prior to re- 
heating at redshift z > 14; and the second generation pop II at z < 10 when 
the virial temperature of the gas clumps reaches 10 4 K and hydrogen line cool- 
ing becomes efficient. Star formation slows in the intermittent epoch due to 
the depletion of H2 by photo-destruction and reheating. In addition, the ob- 
jects which formed pop III stars also initiate pop II sequences when their virial 
temperatures reach 10 4 K through continued mass accretion. 

In resolving the details of a single star forming region in a CDM Universe, 
Abel et al. 2 , 0] implemented a non-equilibrium radiative cooling and chem- 
istry model [1, n9[ together with the hydrodynamics and dark matter equations, 
evolving nine separate atomic and molecular species (H, H + , He, He + , He ++ , 
H , Hj, H2, and e _ ) on nested and adaptively refined numerical grids. They 
follow the collapse and fragmentation of primordial clouds over many decades 
in mass and spatial dynamical range, finding a core of mass ~ 200 Mq forms 
from a halo of about ~ 10 Mq (where a significant number fraction of hydro- 
gen molecules are created) after less than one percent of the halo gas cools by 
molecular line emission. Bromm et al. J43| use a different Smoothed Particle 
Hydrodynamics (SPH) technique and a six species model (H, H + , H~, H^", H2, 
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and e - ) to investigate the initial mass function of the first generation pop III 
stars. They evolve an isolated 3a peak of mass 2 x 1O 6 M0 which collapses at 
rcdshift z ~ 30 and forms clumps of mass 10 2 — 10 3 Afo which then grow by 
accretion and merging, suggesting that the very first stars were massive and in 
agreement with || . 



4.4 Lyman-alpha forest 

The Lyman-alpha forest represents the optically thin (at the Lyman edge) com- 
ponent of Quasar Absorption Systems (QAS), a collection of absorption fea- 
tures in quasar spectra extending back to high redshifts. QAS are effective 
probes of the matter distribution and the physical state of the Universe at early 
epochs when structures such as galaxies are still forming and evolving. Although 
stringent observational constraints have been placed on competing cosmological 
models at large scales by the COBE satellite and over the smaller scales of our 
local Universe by observations of galaxies and clusters, there remains sufficient 
flexibility in the cosmological parameters that no single model has been estab- 
lished conclusively. The relative lack of constraining observational data at the 
intermediate to high redshifts (0 < z < 5), where differences between competing 
cosmological models are more pronounced, suggests that QAS can potentially 
yield valuable and discriminating observational data. 

Several combined N-body and hydrodynamic numerical simulations of the 



Lyman forest have been performed recently ( pi] , |94| , 129 1, for example), and all 
have been able to fit the observations remarkably well, including the column 
density and Doppler width distributions, the size of absorbers |5^], and the line 
number evolution. Despite the fact that the cosmological models and parameters 
are different in each case, the simulations give roughly similar results provided 
that the proper ionization bias is used (6i on = (^bh 2 ) 2 /L, where J7b is the bary- 
onic density parameter, h is the Hubble parameter and T is the photoionization 
rate at the hydrogen Lyman edge). However, see jl5| for a discussion of the 
sensitivity of statistical properties on numerical resolution, and p3| for a sys- 
tematic comparison of five different cosmological models to determine which 
attributes are sensitive physical probes or discriminators of models. A theoret- 
ical paradigm has thus emerged from these calculations in which Lyman-alpha 
absorption lines originate from the relatively small scale structure in pregalac- 
tic or intergalactic gas through the bottom-up hierarchical formation picture in 
CDM-like Universes. The absorption features originate in structures exhibiting 
a variety of morphologies commonly found in numerical simulations (see Fig- 
ure^), including fluctuations in underdense regions, spheroidal minihalos, and 
filaments extending over scales of a few megaparsecs. 
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Figure 6: Distribution of the gas density at redshift z — 3 from a numeri- 
cal hydrodynamics simulation of the Lyman-alpha forest with a CDM spectrum 
normalized to second year COBE observations, Hubble parameter of h = 0.5, a 
comoving box size of 9.6 Mpc, and baryonic density of flh = 0.06 composed of 
76% hydrogen and 24% helium. The region shown is 2.4 Mpc (proper) on a side. 
The isosurfaces represent baryons at ten times the mean density and are color 
coded to the gas temperature ( dark blue = 3 x 10 4 K, light blue = 3x 10 5 K). The 
higher density contours trace out isolated spherical structures typically found at 
the intersections of the filaments. A single random slice through the cube is also 
shown, with the baryonic overdensity represented by a rainbow-like color map 
changing from black (minimum) to red (maximum). The He + mass fraction is 
shown with a wire mesh in this same slice. To emphasize fine structure in the 
minivoids, the mass fraction in the overdense regions has been rescaled by the 
gas overdensity wherever it exceeds unity. 



4.5 Galaxy clusters 

Clusters of galaxies are the largest gravitationally bound systems known to be in 
quasi-equilibrium. This allows for reliable estimates to be made of their mass as 
well as their dynamical and thermal attributes. The richest clusters, arising from 
3cr density fluctuations, can be as massive as 10 14 -10 15 solar masses, and the en- 
vironment in these structures is composed of shock heated gas with temperatures 
of order 10 7 -10 8 degrees Kelvin which emits thermal bremsstrahlung and line 
radiation at X-ray energies. Also, because of their spatial size of ~ 1 h^ 1 Mpc 
and separations of order 50 h^ 1 Mpc, they provide a measure of nonlinearity on 
scales close to the perturbation normalization scale 8/i _1 Mpc. Observations of 
the substructure, distribution, luminosity, and evolution of galaxy clusters are 
therefore likely to provide signatures of the underlying cosmology of our Uni- 
verse, and can be used as cosmological probes in the observable redshift range 
< z < 1. 



4.5.1 Internal structure 



Thomas et al. [ 120 1 investigated the internal structure of galaxy clusters formed 
in high resolution N-body simulations of four different cosmological models, 
including standard, open, and flat but low density Universes. They find that the 
structure of relaxed clusters is similar in the critical and low density Universes, 
although the critical density models contain relatively more disordered clusters 
due to the freeze-out of fluctuations in open Universes at late times. The profiles 
of relaxed clusters are very similar in the different simulations since most clusters 
are in a quasi-equilibrium state inside the virial radius and generally follow 



the universal density profile of Navarro et al. [ 100 1 . There does not appear to 



be a strong cosmological dependence in the profiles as suggested by previous 
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studies of clusters formed from pure power law initial density fluctuations j59| . 
However, because more young and dynamically evolving clusters are found in 
critical density Universes, Thomas et al. suggest that it may be possible to 
discriminate among the density parameters by looking for multiple cores in the 
substructure of the dynamic cluster population. They note that a statistical 
population of 20 clusters could distinguish between open and critically closed 
Universes. 



4.5.2 Number density evolution 

The evolution of the number density of rich clusters of galaxies can be used to 
compute fio and as (the power spectrum normalization on scales of 8 Mpc) 
when numerical simulation results are combined with the constraint ergOf) ~ 
0.5, derived from observed present-day abundances of rich clusters. Bahcall et 
al. [^2) computed the evolution of the cluster mass function in five different 
cosmological model simulations and find that the number of high mass (Coma- 
like) clusters in flat, low as models (i.e., the standard CDM model with cr 8 ~ 0.5) 
decreases dramatically by a factor of approximately 10 3 from z — to z w 0.5. 
For low fio, high as models, the data result in a much slower decrease in the 
number density of clusters over the same redshift interval. Comparing these 
results to observations of rich clusters in the real Universe, which indicate only 
a slight evolution of cluster abundances to redshifts z ~0.5-l, they conclude 
that critically closed standard CDM and Mixed Dark Matter (MDM) models 
are not consistent with the observed data. The models which best fit the data 
are the open models with low bias (f2o = 0.3 ± 0.1 and as = 0.85 ± 0.5), and 
flat low density models with a cosmological constant (fio = 0.34 ± 0.13 and 
O +A=l). 



4.5.3 X-ray luminosity function 

The evolution of the X-ray luminosity function, as well as the number, size and 
temperature distribution of galaxy clusters are all potentially important discrim- 
inants of cosmological models and the underlying initial density power spectrum 
that gives rise to these structures. Because the X-ray luminosity (principally 
due to thermal bremsstrahlung emission from electron/ion interactions in the 
hot fully ionized cluster medium) is proportional to the square of the gas den- 
sity, the contrast between cluster and background intensities is large enough to 
provide a window of observations that is especially sensitive to cluster structure. 
Comparisons of simulated and observed X-ray functions may be used to deduce 
the amplitude and shape of the fluctuation spectrum, the mean density of the 
Universe, the mass fraction of baryons, the structure formation model, and the 
background cosmological model. 

Several groups H, Q have examined the properties of X-ray clusters in 
high resolution numerical simulations of a standard CDM model normalized to 
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COBE. Although the results are very sensitive to grid resolution (see [|15| for 
a discussion of the effects from resolution constraints on the properties of rich 
clusters) , their primary conclusion, that the standard CDM model predicts too 
many bright X-ray emitting clusters and too much integrated X-ray intensity, 
is robust since an increase in resolution will only exaggerate these problems. 
On the other hand, similar calculations with different cosmological models |30|, 
[irif suggest reasonable agreement of observed data with Cold Dark Matter + 
cosmological constant (ACDM), Cold + Hot Dark Matter (CHDM), and Open 
or low density CDM (OCDM) evolutions due to different universal expansions 
and density power spectra. 

4.5.4 SZ effect 

The Sunyaev-Zel'dovich (SZ) effect is the change in energy that CMB photons 
undergo when they scatter in hot gas typically found in cores of galaxy clusters. 
There are two main effects: thermal and kinetic. Thermal SZ is the dominant 
mechanism which arises from thermal motion of gas in the temperature range 
10 7 -10 8 K, and is described by the Compton y parameter 

f n c k B T e 

y = (JT o-dl, (9) 



where <tt = 6.65 x 10 -25 cm 2 is the Thomson cross-section, m c , n c and T e are 
the electron rest mass, density and temperature, c is the speed of light, fee is 
Boltzmann's constant, and the integration is performed over the photon path. 
Photon temperature anisotropies are related to the y parameter by AT/T « — 2y 
in the Rayleigh- Jeans limit. The kinetic SZ effect is a less influential Doppler 
shift resulting from the bulk motion of ionized gas relative to the rest frame of 
the CMB. 



Springel et al. [117] used a Tree/SPH code to study the SZ effects in a 
CDM cosmology with a cosmological constant. They find a mean amplitude for 
thermal SZ (y = 3.8 x 10 -6 ) just below current observed upper limits, and a 



kinetic SZ about 30 times smaller in power. Da Silva et al. |60 compared thermal 
SZ maps in three different cosmologies (low density + A, critical density, and low 
density open model). Their results are also below current limits: j/w4x 10~ 6 
for low density models with contributions from over a broad redshift range z < 5, 
and y ss 1 x 10 -6 for the critical density model with contributions mostly from 
z < 1. However, further simulations are needed to explore the dependence of 
the SZ effect on microphysics, i.e., cooling, star formation, supernovae feedback. 

4.6 Cosmological sheets 

Cosmological sheets, or pancakes, form as overdense regions collapse prefer- 



entially along one axis. Originally studied by Zel'dovich [ 128 1 in the context 
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of neutrino-dominated cosmologies, sheets are ubiquitous features in nonlin- 
ear structure formation simulations of CDM-like models with baryonic fluid, 
and manifest on a spectrum of length scales and formation epochs. Gas col- 
lapses gravitationally into flattened sheet structures, forming two plane parallel 
shock fronts that propagate in opposite directions, heating the infalling gas. The 
heated gas between the shocks then cools radiatively and condenses into galactic 
structures. Sheets are characterized by essentially five distinct components: the 
preshock inflow, the postshock heated gas, the strongly cooling/recombination 
front separating the hot gas from the cold, the cooled postshocked gas, and the 
unshocked adiabatically compressed gas at the center. Several numerical calcu- 
2p| have been performed of these systems which include baryonic 
fluid with hydrodynamical shock heating, ionization, recombination, dark mat- 
ter, thermal conductivity, and radiative cooling (Compton, bremsstrahlung, and 
atomic line cooling), in both one and two spatial dimensions to assert the sig- 
nificance of each physical process and to compute the fragmentation scale. See 
also ]I3 where fully general relativistic numerical calculations of cosmological 
sheets are presented in plane symmetry, including relativistic hydrodynamical 
shock heating and consistent coupling to spacetime curvature. 



lations M2, 113, 



Figure 7: Two different model simulations of cosmological sheets are presented: 
a six species model including only atomic line cooling (left), and a nine species 
model including also hydrogen molecules (right). The evolution sequences in the 
images show the baryonic overdensity and gas temperature at three redshifts fol- 
lowing the initial collapse at z — 5. In each figure, the vertical axis is 32 kpc 
long (parallel to the plane of collapse) and the horizontal axis extends to 4 Mpc 
on a logarithmic scale to emphasize the central structures. Differences in the 
two cases are observed in the cold pancake layer and the cooling flows between 
the shock front and the cold central layer. When the central layer fragments, the 
thickness of the cold gas layer in the six (nine) species case grows to 3 (0.3) kpc 
and the surface density evolves with a dominant transverse mode corresponding 
to a scale of approximately 8 (1) kpc. Assuming a symmetric distribution of 
matter along the second transverse direction, the fragment masses are approxi- 
mately 10 7 fl0 5 J solar masses. 

In addition, it is well known that gas which cools to 1 eV through hydrogen 
line cooling will likely cool faster than it can recombine. This nonequilibrium 
cooling increases the number of electrons and ions (compared to the equilibrium 
case) which, in turn, increases the concentrations of H~ and , the intermedi- 
aries that produce hydrogen molecules H 2 . If large concentrations of molecules 
form, excitations of the vibrational/rotational modes of the molecules can ef- 
ficiently cool the gas to well below 1 eV, the minimum temperature expected 
from atomic hydrogen line cooling. Because the gas cools isobarically, the re- 
duction in temperature results in an even greater reduction in the Jeans mass, 
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and the bound objects which form from the fragmentation of H2 cooled cosmo- 
logical sheets may be associated with massive stars or star clusters. Anninos 
and Norman Jl6| have carried out ID and 2D high resolution numerical calcu- 
lations to investigate the role of hydrogen molecules in the cooling instability 
and fragmentation of cosmological sheets, considering the collapse of perturba- 
tion wavelengths from 1 Mpc to 10 Mpc. They find that for the more ener- 
getic (long wavelength) cases, the mass fraction of hydrogen molecules reaches 
«H2/^h ~3x 10~ 3 , which cools the gas to 4 x 10~ 3 eV and results in a frag- 
mentation scale of 9 x 10 3 M Q . This represents reductions of 50 and 10 3 in 
temperature and Jeans mass respectively when compared, as in Figure to the 
equivalent case in which hydrogen molecules were neglected. 

However, the above calculations neglected important interactions arising 
from self-consistent treatments of radiation fields with ionizing and photo-disso- 



ciating photons and self-shielding effects. Susa and Umemura [ 1 18 1 studied the 



thermal history and hydrodynamical collapse of pancakes in a UV background 
radiation field. They solve the radiative transfer of photons together with the 
hydrodynamics and chemistry of atomic and molecular hydrogen species. Al- 
though their simulations were restricted to one-dimensional plane parallel sym- 
metry, they suggest a classification scheme distinguishing different dynamical 
behavior and galaxy formation scenarios based on the UV background radia- 
tion level and a critical mass corresponding to 1 — 2a density fluctuations in a 
standard CDM cosmology. These level parameters distinguish galaxy formation 
scenarios as they determine the local thermodynamics, the rate of H2 line emis- 
sions and cooling, the amount of starburst activity, and the rate and mechanism 
of cloud collapse. 
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5 Conclusion 



This review is intended to provide a flavor of the variety of numerical cos- 
mological calculations performed of different phenomena occurring throughout 
the history of our Universe. The topics discussed range from the strong field 
dynamical behavior of spacetime geometry at early times near the Big Bang 
singularity and the epoch of inflation, to the late time evolution of large scale 
matter fluctuations and the formation of clusters of galaxies. Although a com- 
plete, self-consistent, and accurate description of our Universe is impractical 
considering the complex multiscale and multiphysics requirements, a number 
of enlightening results have been demonstrated through computations. For ex- 
ample, both monotonic AVTD and chaotic oscillatory BLK behavior have been 
found in the asymptotic approach to the initial singularity in a small set of 
inhomogencous Bianchi and Gowdy models, though it remains to be seen what 
the generic behavior might be in more general multidimensional spacetimes. 
Numerical calculations suggest that scalar fields play an important complicated 
role in the nonlinear or chaotic evolution of cosmological models with conse- 
quences for the triggering (or not) of inflation and the subsequent dynamics 
of structure formation. It is possible, for example, that these fields can influ- 
ence the details of inflation and have observable ramifications as fractal patterns 
in the density spectrum, gravitational waves, galaxy distribution, and cosmic 
microwave background anisotropies. All these effects require further studies. 
Numerical simulations have been used to place limits on curvature anisotropics 
and cosmological parameters at early times by considering primordial nucle- 
osynthesis in anisotropic and inhomogeneous cosmologies. Finally, the large 
collection of calculations performed of the post-recombination epoch (for ex- 
ample, cosmic microwave, gravitational lensing, Lyman-alpha absorption, and 
galaxy cluster simulations) have placed strong constraints on the standard model 
parameters and structure formation scenarios when compared to observations. 
Considering the range of models consistent with inflation, the preponderance 
of observational, theoretical and computational data suggest a best fit model 
that is spatially flat with a cosmological constant and a small tilt in the power 
spectrum. 

Obviously many fundamental issues remain unresolved, including the back- 
ground or topology of the cosmological model which best describes our Universe, 
the generic singularity behavior, the dynamics of inflaton fields, the imprint 
of complex interacting scalar fields, the fundamental nonlinear curvature and 
gravitational wave interactions, the correct structure formation scenario, and 
the origin and spectrum of primordial fluctuations, for example, are uncertain. 
However the field of numerical cosmology has matured in the development of 
computational techniques, the modeling of microphysics, and in taking advan- 
tage of current computing technologies, to the point that it is now possible to 
perform high resolution multiphysics simulations and reliable comparisons of 
numerical models with observed data. 
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6 Appendix: Basic Equations and Numerical 
Methods 



Some basic equations relevant for fully relativistic as well as perturbative cosmo- 
logical calculations are summarized in this section, including the complete Ein- 
stein equations, choices of kinematical conditions, initial data constraints, stress- 
energy-momentum tensors, dynamical equations for various matter sources, and 
the Newtonian counterparts on background isotropic models. References to nu- 
merical methods are also provided. 

6.1 The Einstein equations 
6.1.1 ADM formalism 

There are many ways to write the Einstein equations. The most common is 
the ADM or 3 + 1 form [glj which decomposes spacetime into layers of three- 
dimensional space- like hypersurfaces, threaded by a time-like normal congruence 
n M = (1, —(3 l )/a. The general spacetime metric is written as 

ds 2 = (-a 2 + fafi^dt 2 + 2/3 i dx i dt + ^jdx l dx\ (10) 

where a and /3 Z are the lapse function and shift vector respectively, and 7^ is the 
spatial 3-metric. The lapse defines the proper time between consecutive layers of 
spatial hypersurfaces, the shift propagates the coordinate system from 3-surface 
to 3-surface, and the induced 3-metric is related to the 4-metric via 7^ — 
9fiv + n^riy. The Einstein equations are written as four constraint equations, 

(3) i? + K 2 - /\ , ,/\ ■' = levrGpH, (11) 
Vi (K ij - f j K) = 8nGs j , (12) 

twelve evolution equations, 

dt-fij = Cpji-j - 2aK^, 

d t Kij = CpKi-j - V;VjcH- q 3 \ 
WRij - 2K lk K* + KK t j - 8irG ( Sij - ^ Sllj + ^h7« 

and four kinematical or coordinate conditions for the lapse function and shift 



vector that can be specified arbitrarily (see § 6.2). Here, 

Cpjij =Vi/3,- +Vj/3j, 
CpKij = (i k ^ k K l3 + K lk S7 J p k + K kj Vif3 k , 

where Vi is the spatial covariant derivative with respect to 7y , *- 3 ^ Rij the spatial 
Ricci tensor, K the trace of the extrinsic curvature Kij, and G is the gravita- 
tional constant. The matter source terms pn, s 3 , Sy and s = s\ as seen by the 
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observers at rest in the time slices are obtained from the appropriate projections 



p H = n^n v T llv , (15) 



•s 



= - 7 fn%„, (16) 
s y = lilj T ^ ( 17 ) 

for the energy density, momentum density and spatial stresses, respectively. 
Here c = 1, and Greek (Latin) indices refer to 4(3)-dimensional quantities. 
It is worth noting that several alternative formulations of Einstein's equa- 



tions have been suggested, including hyperbolic systems [ 105 1 which have nice 
mathematical properties, and conformal traceless systems [115, which make 
use of a conformal decomposition of the 3-metric and trace-free part of the ex- 
trinsic curvature. Introducing 7^ = e~ 4 ^7y with e 4 ^ = 7 1 / 3 so that the deter- 
minant of jij is unity, and Aij — e~^Aij, evolution equations can be written in 
the conformal traceless system for -0, 7^, K, Aij and the conformal connection 
functions, though not all of these variables are independent. However, it is not 
yet entirely clear which of these methods is best suited for generic problems. For 
example, hyperbolic forms are easier to characterize mathematically than ADM 
and may potentially be more stable, but can suffer from greater inaccuracies by 
introducing additional equations or higher order derivatives. Conformal treat- 
ments are considered to be generally more stable J 28J , but can be less accurate 
than traditional ADM for short term evolutions |5]. 

Many numerical methods have been used to solve the Einstein equations, 
including variants of the leapfrog scheme, the method of McCormack, the two- 
step Lax-Wendroff method, and the iterative Crank-Nicholson scheme, among 
others. For a discussion and comparison of the different methods, the reader 
is referred to |39f| , where a systematic study was carried out on spherically 
symmetric black hole spacetimes using traditional ADM, and to j|§ |, [TJ (and 
references therein) which discuss the stability and accuracy of hyperbolic and 
conformal treatments. 



6.1.2 Symplectic formalism 

A different approach to conventional (i.e., 3 + 1 ADM) techniques in numerical 
cosmology has been developed by Berger and Moncrief J38|. For example, they 
consider Gowdy cosmologies on T 3 x R with the metric 

ds 2 = e - A /V /2 (~e- 2T dr 2 + de 2 ) + 

e~ T [e p da 2 + 2e p Qdad5 + (e p Q 2 + e~ p ) dS 2 ] , (18) 

where A, P and Q are functions of 9 and r, and the coordinates are bounded by 
< (9,(7,5) < 2ir. The singularity corresponds to the limit t — > 00. For small 
amplitudes, P and Q may be identified with + and x polarized gravitational 
wave components and A with the background cosmology through which they 
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propagate. An advantage of this formalism is that the initial value problem 
becomes trivial since P, Q and their first derivatives may be specified arbitrarily 
(although it is not quite so trivial in more general spacetimes). 

Although the resulting Einstein equations can be solved in the usual space- 
time discretization fashion, an interesting alternative method of solution is the 
symplectic operator splitting formulation []38| , |9lj| founded on recognizing that 
the second order equations can be obtained from the variation of a Hamiltonian 
decomposed into kinetic and potential subhamiltonians, 

H = Hx+H 2 = ^i d6(ir% + e~ 2P ir 2 Q ) + dOe^ (P 2 e + e 2P Q 2 g ) . 

(19) 

The symplectic method approximates the evolution operator by 

e HAr = e ff 2 Ar/2 e ff I Ar e ff 2 Ar/2 + 0( At )3 ; (2Q) 

although higher order representations are possible. If the two Hamiltonian com- 
ponents H\ and H 2 are each integrable, their solutions can be substituted di- 
rectly into the numerical evolution to provide potentially more accurate solu- 
tions with fewer time steps f^f . This approach is well-suited for studies of 
singularities if the asymptotic behavior is determined primarily by the kinetic 
subhamiltonian, a behavior referred to as Asymptotically Velocity Term Domi- 



nated (see § |3.1.2| and |3J). 

Symplectic integration methods are applicable to other spacetimes. For ex- 
ample, Berger et al. |}5| developed a variation of this approach to explicitly take 
advantage of exact solutions for scattering between Kasner epochs in Mixmaster 
models. Their algorithm evolves Mixmaster spacetimes more accurately with 
larger time steps than previous methods. 

6.2 Kinematic conditions 

For cosmological simulations, one typically takes the shift vector to be zero, 
hence £/37y = CpKij = 0. However, the shift can be used advantageously 
in deriving conditions to maintain the 3-metric in a particular form, and to 



simplify the resulting differential equations |54|, pa]. See also reference [114 
describing an approximate minimum distortion gauge condition used to help 
stabilize simulations of general relativistic binary clusters and neutron stars. 

Several options can be implemented for the lapse function, including geodesic 
(a = l), algebraic, and mean curvature slicing. The algebraic condition takes 
the form 

a = F 1 {x»)F 2 (j), (21) 

where F\{x^) is an arbitrary function of coordinates 2;^, and ^2(7) is a dy- 
namic function of the determinant of the 3-metric. This choice is computation- 
ally cheap, simple to implement, and certain choices of F 2 (i.e., 1 + In 7) can 
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mimic maximal slicing in its singularity avoidance properties [[7| . On the other 
hand, numerical solutions derived from harmonically-sliced foliations can exhibit 
pathological gauge behavior in the form of coordinate "shocks" or singularities 
which will affect the accuracy, convergence and stability of solutions || |74| . 
Also, evolutions in which the lapse function is fixed by some analytically pre- 
scribed method (either geodesic or near-geodesic) can be unstable, especially 
for sub- horizon scale perturbations || . 

The mean curvature slicing equation is derived by taking the trace of the 
extrinsic curvature evolution equation Jl3|), 

V J V. ( a = a [KijK ij + 4ttG (p H + *)] + (3 l V t K - d t K, (22) 

and assuming K = K(t), which can either be specified arbitrarily or determined 
by imposing a boundary condition on the lapse function after solving (|22| ) for the 
quantity a/d t K |5j|. It is also useful to consider replacing d t K in equation ( f22|) 
with an exponentially driven form as suggested by Eppley Q , to reduce gauge 
drifting and numerical errors in maximal and mean curvature Q sliced 
spacetimes. The mean curvature slicing condition is the most natural one for 
cosmology as it foliates homogeneous cosmological spacetimes with surfaces of 
homogeneity. Also, since K represents the convergence of coordinate curves 
from one slice to the next and if it is constant, then localized caustics (crossing 
of coordinate curves) and true curvature singularities can be avoided. Whether 
general inhomogeneous spacetimes can be foliated with constant mean curvature 
surfaces remains unknown. However, for Gowdy spacetimes with two Killing 
fields and topology T 3 x R, Isenberg and Moncrief Q proved that such foliations 
do exist and cover the entire spacetime. 



6.3 Sources of matter 
6.3.1 Cosmological constant 

A cosmological constant is implemented in the 3 + 1 framework simply by intro- 
ducing the quantity — A/(8itG) as an effective isotropic pressure in the stress- 
energy tensor 

T u.v = -g^£ 6V- ( 23 ) 
The matter source terms can then be written 



A 

Z 
A 

with s l = 0. 



Hi, (25) 
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6.3.2 Scalar field 

The dynamics of scalar fields is governed by the Lagrangian density 

C = ~V=9 W lv ^ v + £Rf{<f>) + W{4>)\ , (26) 

where R is the scalar Riemann curvature, V((f>) is the interaction potential, 
f{4>) is typically assumed to be /(</>) = </> 2 , and £ is the field-curvature coupling 
constant (£ = for minimally coupled fields and £ = 1/6 for conformally coupled 
fields). Varying the action yields the Klein-Gordon equation 

g^k^-tRt- W£) = 0, (27) 

for the scalar field and 

T MV = (1 - 20^.„ + fa - 

-2e# i|iV + Wg^hrx + £G^ 2 - g^V(<j>), (28) 

for the stress-energy tensor, where = i? M „ — g UJJ R/2. 
For a massive, minimally coupled scalar field pdj j 

= 0;^</> ;i / - ^g^g^^-.p^-a - gn V V(4>), (29) 



and 



Ph = ^7 ij >;^d + ^ 2 + V(<f>), (30) 

Sj = -T)(j>-i, (31) 

*« = 7« (-^7 fe V;fc0;i + ^ 2 - + (32) 



where 



r? = n^= i(9 t -/3 fc a fe )0, (33) 

n M = (1,— and V r (<^) is a general potential which, for example, can be 

set to V = A</> 4 in the chaotic inflation model. The covariant form of the scalar 
field equation ( p7| ) can be expanded as in to yield 

- (ft - P k d k ) V = ^diiy/rf'drf) + --.'u),,u),o + K V - B^V((f>), (34) 
a a 

which, when coupled to (p3|), determines the evolution of the scalar field. 
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6.3.3 Collisionless dust 

The stress-energy tensor for a fluid composed of collisionless particles (or dark 
matter) can be written simply as the sum of the stress-energy tensors for each 
particle [p5| , 

Tfiu = mnu^u u , (35) 

where m is the rest mass of the particles, n is the number density in the comoving 
frame, and u M is the 4-velocity of each particle. The matter source terms are 

Ph = ^mn(aii ) 2 , (36) 
8} = 'S~^mnUi(au°), (37) 
Sij = mnUiUj . (38) 

There are two conservation laws: the conservation of particles V ^(nu^) = 0, and 
the conservation of energy-momentum V M T M1/ = 0, where V M is the covariant 
derivative in the full 4-dimcnsional spacetime. Together these conservation laws 
lead to it^V^u^ = 0, the geodesic equations of motion for the particles, which 
can be written out more explicitly in the computationally convenient form 

dx l q la u n , . 

*" = V' < 39 > 

duj __ u a U0djg a P 

dt ~ 2u° ' ( ' 

where x 1 is the coordinate position of each particle, u° is determined by the 
normalization u^u^ = — 1, 

± = v*d li = a t + v i a i (4i) 

is the Lagrangian derivative, and = u^/u is the "transport" velocity of the 
particles as measured by observers at rest with respect to the coordinate grid. 

6.3.4 Ideal gas 

The stress-energy tensor for a perfect fluid is 

Tpv = phu^Uv + Pg^ v , (42) 

where g^ v is the 4-metric, h = 1 + e + Pjp is the relativistic enthalpy, and e, P, 
p and iiy. are the specific internal energy (per unit mass), pressure, rest mass 
density and four-velocity of the fluid. Defining 



— 1/2 

*" = ou° = (1 + u'ui) 1/2 = ( 1 - ™P\ , (43) 
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as the generalization of the special relativistic boost factor, the matter source 
terms become 



Si 
Sij 



phu 2 - P, 
phuui, 

Pltj + phuiUj. 



(44) 
(45) 
(46) 



The hydrodynamics equations are derived from the normalization of the 4- 



velocity, u^u^ 



T, the conservation of baryon number, V At (pw' i ) = 0, and 



the conservation of energy-momentum, V U T IJ ' V = 0. The resulting equations 



can be written in flux conservative form as [126 

dD d(DV 

d(EV l ) 



dE 
~dt 

dt 



dx l 
d{StV*) 



dx % 

~dt 

S»S" d 9llu 



0. 



P- 



P 



d(W V 1 ) 
dx i 



.dP 



= 0, 
V 



(47) 
(48) 
(49) 



"dx* 

Wphu i: V 1 — u l /u°, and g 
J—g = a-y/j. A prescription for 
1)E/W = (T - l)pe for an ideal 



dxi 2S° dx l 
where W = >J=gu , D = Wp, E — Wpe, 5, 
is the determinant of the 4-metric satisfying ^ 
specifying an equation of state (e.g., P = (T 
gas) completes the above equations. 

When solving equations ([l7|, |4^, ^) , an artificial visc osity method is needed 
to handle the formation and propagation of shock fronts [126, 72, 7j|. Although 
these methods are simple to implement and are computationally efficient, they 
are inaccurate for ultrarelativistic flows with very high Lorentz factors. On 
the other hand, a number of different formulations ]67|] of these equations have 
been developed to take advantage of the hyperbolic and conservative nature 
of the equations in using high resolution shock capturing schemes (although 
strict conservation is only possible in flat spacetimes - curved spacetimes exhibit 
source terms due to geometry). Such high resolution Godunov techniques [107, 
pi| provide more accurate and stable treatments in the highly relativistic regime. 
A particular formulation due to [M is the following: 



<Vr U{w) dy/=g F l (w) 



dt 



dx l 



= a/ 3 ;? S(w), 



where 



S(w) 
F l (w) 



Dlv 1 - 1 - \SAv l -^) + P8),{E-D)[v*-^- 



a 



Pv l 



(50) 

(51) 
(52) 



and w — (p, vi , e) 
D = pW, v' 



U(w) = (D,Si,E- D), E = phW 2 - P, Sj = phW 2 
u l /{au°) + j3 l /a, and W = au° = (1 - 7y<uV)- 1 / 2 . 
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6.4 Constrained nonlinear initial data 



One cannot take arbitrary data to initiate an evolution of the Einstein equa- 
tions. The data must satisfy the constraint equations ( |li"| ) and York [ 127 1 



developed a procedure to generate proper initial data by introducing conformal 
transformations of the 3- metric 7y = ip^^ij, the trace-free momentum compo- 
nents = K % i — 7 1 - 7 K/3 = i(;~ 1Q A l: > , and matter source terms s l = ?/; _10 s 4 
and ph = "0 ™/5h, where n > 5 for uniqueness of solutions to the elliptic equa- 
tion ( |53| ) below. In this procedure, the conformal (or "hatted") variables are 
freely specifiable. Further decomposing the free momentum variables into trans- 
verse and longitudinal components A^ — A % i + (/«;)*■? , the Hamiltonian and 
momentum constraints are written as 



ViV*V - + IkjA^^- 7 - + 2nGpiP 5 - n = 0, (53) 

o o 12 



(V,-V J u;) 1 + (V 3 w 3 \ + mw j - \^ l K - 8nGs % = 0, (54) 
3 V / J 3 

where the longitudinal part of A 13 is reconstructed from the solutions by 

{Iwfi = vV + vV - \f 3 V k w k . (55) 

The transverse part of A^ is constrained to satisfy VjAl J = Ajj = 0. 

Equations (|5^) and (|54|) form a coupled nonlinear set of elliptic equations 
which must be solved iteratively, in general. The two equations can, however, 
be decoupled if a mean curvature slicing (K = K(t)) is assumed. Given the 
free data K, %j, s l and p, the constraints are solved for A l £ , (iw) 1 ^ and ip. 
The actual metric 7^- and curvature ify are then reconstructed by the corre- 
sponding conformal transformations to provide the complete initial data. Refer- 
ence H describes a procedure using York's formalism to construct parametrized 
inhomogeneous initial data in freely specifiable background spacetimes with 
matter sources. The procedure is general enough to allow gravitational wave 
and Coulomb nonlinearities in the metric, longitudinal momentum fluctuations, 
isotropic and anisotropic background spacetimes, and can accommodate the 
conformal-Newtonian gauge to set up gauge invariant cosmological perturba- 
tion solutions as free data. 



6.5 Newtonian limit 

The Newtonian limit is defined by spatial scales much smaller than the horizon 
radius, peculiar velocities small compared to the speed of light, and a gravita- 
tional potential that is both much smaller than unity (in geometric units) and 
slowly varying in time. A comprehensive review of the theory of cosmological 
perturbations can be found in |99fl . 
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6.5.1 Dark and baryonic matter equations 

The appropriate perturbation equations in this limit are easily derived for a 
background FLRW expanding model, assuming a metric of the form 

ds 2 = -(1 + 2$)dt 2 + a(t) 2 (l - 2$)j ij dx i dx j , (56) 

where 

™ = 3 ( x + ^r) - 2 ( 57 ) 

and fc = — 1, 0, +1 for open, flat and closed Universes. Also, a = 1/(1 + 2;) is the 
cosmological scale factor, z is the redshift, and 4> is the comoving inhomogeneous 
gravitational potential. 

The governing equations in the Newtonian limit are the hydrodynamic con- 
servation equations, 

f + + sifc-0, (58) 

the geodesic equations for collisionless dust or dark matter (in comoving coor- 
dinates), 

(61) 

c(t a a z ox 1 

Poisson's equation for the gravitational potential, 

V 2 $- 4irGa 2 (p h + p d -p Q ), (63) 

and the Friedman equation for the cosmological scale factor, 



da 

dt ~ Ho 



O m (--l) + fi A (a 2 -l) + l 



1/2 



(64) 



Here pd, pb, P and e are the dark matter density, baryonic density, pressure 
and internal energy density in the proper reference frame, x % and v h are the 
baryonic comoving coordinates and peculiar velocities, x d and v d are the dark 
matter comoving coordinates and peculiar velocities, po = 311^0 /(8nGa 3 ) is 
the proper background density of the Universe, O is the total density parame- 
ter, O m = fl b + fid is the density parameter including both baryonic and dark 
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matter contributions, fl\ = A/(3Hq) is the density parameter attributed to 
the cosmological constant A, Ho = 100 h km s" 1 Mpc -1 is the present Hub- 
ble constant with 0.5 < h < 1, and Scooi represents microphysical radiative 
cooling and heating rates which can include Compton cooling (or heating) due 
to interactions of free electrons with the CMBR, bremsstrahlung, and atomic 
and molecular line cooling. Notice that 'tilded' ('non-tilded') variables refer to 
proper (comoving) reference frame attributes. 

An alternative total energy conservative form of the hydrodynamics equa- 
tions that allows high resolution Godunov-type shock capturing techniques is 

d{p h v 3 h ) 13 ~i ~j . * n . a , pi, d& , , 

—QT- + ad^ iPhV * < +p6i ^ + a Ph < + TfeJ = °' (66) 

M + aW^ +jn ^ + ~a ad^ = (67) 

with the corresponding particle and gravity equations 

# - !• <«> 
§ - 4« - \t < 69 » 

V 2 $ = — Ob + Pd-po), (70) 

where pb is the comoving density, p — a3 Po, w b is the proper frame peculiar 
velocity, p is the comoving pressure, E = phv\ /2 +p/ (7— 1) is the total peculiar 
energy per comoving volume, and $ is the gravitational potential. 

These equations are easily extended JTsf] to include reactive chemistry of nine 
separate atomic and molecular species (H, H + , He, He + , He ++ , H~, Hj, H 2 , and 
e _ ), assuming a common flow field, supplementing the total mass conservation 
equation (|5^) with 



dt dx % 



a 

i I 



for each of the species, and including the effects of non-equilibrium radiative 
cooling and consistent coupling to the hydrodynamics equations. The ku(T) 
are rate coefficients for the two body reactions and are tabulated functions of 
the gas temperature T. The Ii are integrals evaluating the photoionization and 
photodissociation of the different species. For a comprehensive discussion of the 
cosmologically important chemical reactions and reaction rates, see reference Q. 
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Many numerical techniques have been developed to solve the hydrodynamic 
and collisionless particle equations. For the hydrodynamic equations, the meth- 
ods range from Lagrangian SPH algorithms p4], [76[ |103| to Eulcrian finite 
difference techniques on static meshes [110, |104| , nested grids |L7|, moving 
meshes JfO] ]) and adaptive mesh refine men t |47j[. For the dark matter equations, 
the canonical choices arc trcccodes [1231 or PM and P 3 M methods [[781 p2|, 



although many variants have been developed to optimize computational per- 
formance and accuracy, including grid and particle refinement methods (see 
references cited in p8f ). An efficient method for solving non-equilibrium, multi- 
species chemical reactive flows together with the hydrodynamic equations in a 
background FLRW model is described in jl], [li| . 



The reader is referred to |83 



for thorough comparisons of different nu- 
merical methods applied to problems of structure formation. Reference Q 
compares (by binning data at different resolutions) the statistical performance 
of five codes (three Eulerian and two SPH) on the problem of an evolving CDM 
Universe on large scales using the same initial data. The results indicate that 
global averages of physical attributes converge in rebinned data, but that some 
uncertainties remain at small levels. |68|] compares twelve Lagrangian and Eule- 
rian hydrodynamics codes to resolve the formation of a single X-ray cluster in a 
CDM Universe. The study finds generally good agreement for both dynamical 
and thermodynamical quantities, but also shows significant differences in the 
X-ray luminosity, a quantity that is especially sensitive to resolution flTa . 



6.5.2 Linear initial data 



The standard Zel'dovich solution [128, 62 can be used to generate initial condi- 
tions satisfying observed or theoretical power spectra of matter density fluctu- 
ations. Comoving physical displacements and velocities of the collisionless dark 
matter particles are set according to the power spectrum realization 



Sp 
P 



(k) 



oc k n T 2 (k), 



(72) 



where the complex phases are chosen from a gaussian random field, T(k) is a 
transfer function [ p5[ appropriate to a particular structure formation scenario 
(e.g., CDM), and n = 1 corresponds to the Harrison-Zel'dovich power spectrum. 
The fluctuations are normalized with top hat smoothing using 

1 r°° 

4irk 2 P{k)W 2 (k)dk, (73) 



9 

^8 



b 2 



where b is the bias factor chosen to match present observations of rms density 
fluctuations in a spherical window of radius Rh = 8 h^ 1 Mpc. Also, P(k) is the 
Fourier transform of the square of the density fluctuations in equation ( [72] ) , and 

W(k) = jj^-^ (sm(kR h ) - (kR h ) cos(fci? h )) (74) 
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is the Fourier transform of a spherical window of radius Rh- 

Overdensity peaks can be filtered on specified spatial or mass scales by Gaus- 
sian smoothing the random density field pa] 



^»-p4p/7 M - p (- ls W l! ) A ' ,75) 

on a comoving scale Rf centered at r = r Q (for example, Rf = 5/i _1 Mpc 
with a filtered mass of Mf ~ 10 15 M Q over cluster scales). Na peaks are gen- 
erated by sampling different random field realizations to satisfy the condition 
v = <r(r )/cr(Rf) = N, where <r(R{) is the rms of Gaussian filtered density 
fluctuations over a spherical volume of radius R{ . 

Bertschinger [[lO] has provided a useful and publicly available package of pro- 
grams called COSM1CS for computing transfer functions, CMB anisotropics, 
and gaussian random initial conditions for numerical structure formation cal- 
culations. The package solves the coupled linearized Einstein, Boltzman, and 
fluid equations for scalar metric perturbations, photons, neutrinos, baryons, and 
collisionless dark matter in a background isotropic Universe. It also generates 
constrained or unconstrained matter distributions over arbitrarily specifiable 
spatial or mass scales. 
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